subroutine initialization
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M
      use cophys_com_M
      use ctemp_com_M
      use blcom_com_M
      use objects_com_M
      use controllo_M
      use celeste3d
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------


     open(unit=5, file='input', position='asis')
!ll   1. open i/o devices:
!        ----------------
!      call link ("unit22 = (input   ,open),
!     .            unit9  = (out1 ,create,text),
!     .            unit6  = (out2 ,create,text),
!     .            unit99 = (film    ,create,text),
!     .            unit15 = (pldata,create,len=10b),
!     .            unit16 = (box     ,create,len=10b),
!     .            //")
!

!
!
!l    krd    ...device number for read input cards
!l    kpt    ...device number for print
!l    kpr    ...device number for print (switched)
!l    kpl    ...device number for storing data on disk
!l    kfm    ...device number for film/microfiche
!l    ktpin  ...device number for binary tape read
!l    ktpout ...device number for binary tape write
      krd = 5
      kpt = 6
      kpl = 15
      kfm = 59
      ktpin = 7
      ktpout = 7
!
!ll   2. initialize run:
!        --------------
!
!     initialize control counters
!
      sok = 0
      sko = 0
      cok = 0
      cko = 0
!

 	  pixxsp2 = 0.0
          pixysp2 = 0.0
          pixzsp2 = 0.0
          piyysp2 = 0.0
          piyzsp2 = 0.0
          pizzsp2 = 0.0
 	  pixxsp1 = 0.0
          pixysp1 = 0.0
          pixzsp1 = 0.0
          piyysp1 = 0.0
          piyzsp1 = 0.0
          pizzsp1 = 0.0
!
          densp1=0d0
		  uxsp1=0d0
		  uysp1=0d0
		  uzsp1=0d0
!
          densp2=0d0
		  uxsp2=0d0
		  uysp2=0d0
		  uzsp2=0d0
!
          exavp = 0.0
          eyavp = 0.0
          ezavp = 0.0
!
          joldx = 0d0
          joldy = 0d0
          joldz = 0d0
!
          avg=0d0
!
          potavp = 0.0
!
!     open graphics library
!
      call gplot ('u', celplt, 6)
      call libdisp

      write (0, *) ' celest3d: starting initiv'

      call initiv

      if (restart) then
         call restartin
         call plotinit (ibp2, jbp2, kbp2, nrg)
         ncyc1 = 1
         ncyc2 = ncyc_restart
         write (*, *) 'ncyc1,2', ncyc1, ncyc2
         call date_and_time(ddate,hour)
         write(*,*)'date (ccyymmdd) = ',ddate
         write(*,*)'time (hhmmss.sss) = ',hour
         return
      endif
      call plotinit (ibp2, jbp2, kbp2, nrg)
!
!     *******************************************************
!
!     calculate the gradient of the magnetic field strength only if
!     the drift approximation is used
!
!     **********************************************************
      if (anydrift) then
!
         do k = 1, kbp2
            do j = 1, jbp2
               i = 1
               if (ibp2 > 0) then
!
                  absb(1+(j-1)*jwid+(k-1)*kwid:(ibp2-1)*iwid+1+(j-1)*jwid+(k-1)&
                     *kwid:iwid) = 0.0
!
                  i = ibp2 + 1
                  ijk = (ibp2 - 1)*iwid + 1 + (j - 1)*jwid + (k - 1)*kwid
               endif
            end do
         end do
!
         write (0, *) 'celest3d: starting absolute'
         call absolute (ncells, ijkcell, bxn, byn, bzn, absb)
!c
         wk = 0.
         absbc = 0.
         call bcphi (ibp1, jbp1, kbp1, iwid, jwid, kwid, wkr, wkr, wkf, wke, &
            absbc, absb, periodicx)
!
!
         call gradf (nvtx, ijkvtx, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, &
            c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, &
            c3z, c4z, c5z, c6z, c7z, c8z, absb, gradxb, gradyb, gradzb)
!
         call axisgrad (ibp1, jbp1, iwid, jwid, kwid, gradxb, gradyb, gradzb)
!
!
         do n = 1, nvtx
!
            ijk = ijkvtx(n)
!
            rvvol = 1./vvol(ijk)
!
            gradxb(ijk) = gradxb(ijk)*rvvol
            gradyb(ijk) = gradyb(ijk)*rvvol
            gradzb(ijk) = gradzb(ijk)*rvvol
!
         end do
!
      endif
!
!
!     ************************************************
!
!      initialize the particles
!
!     ************************************************
!
      write (0, *) ' celest3d: starting parset'
      call current
      call parset
      if (nrg_obj > 0) call parset_obj
!
      where = 0.
!      call counter(ncells,ijkcell,iphead,link,number,npart,where)
!
!     **************************************************
!
!     interpolate particle quantities to the grid
!
!     **************************************************
!
      write (0, *) ' celest3d: starting parcelc'
      call parcelc
!
!
      transpos = 1.
!
      call parcelv (ncells, ijkcell, nvtx, ijkvtx, nsp, iwid, jwid, kwid, ibar&
         , jbar, kbar, dt, ijkctmp, anydrift, eandm, itdim, iphead, iphd2, link&
         , transpos, qom, ico, qpar, pxi, peta, pzta, up, vp, wp, pmu, bxv, byv&
         , bzv, wate, uptilde, vptilde, wptilde, mask, pixx, pixy, pixz, piyy, &
         piyz, pizz, qdnv, number, j12x, j12y, j12z, jmag, jxtilde, jytilde, &
         jztilde, jxs, jys, jzs,  pixysp2, piyzsp2, pixzsp2, pixxsp2, piyysp2, pizzsp2, &
         pixysp1, piyzsp1, pixzsp1, pixxsp1, piyysp1, pizzsp1, uxsp1,uysp1,uzsp1,densp1,uxsp2,uysp2,uzsp2,densp2, 0d0)
!
!
!     ******************************************************
!
!     interpolate object quantities to the grid
!
!     *****************************************************
!
      if (nrg_obj > 0) then
!
         write (*, *) 'celest3d: starting parcelv_obj'
!
         call parcelv_obj (ncells, ijkcell, nvtx, ijkvtx, nrg_obj, iwid, jwid, &
            kwid, ijkctmp, itdim, iphead_obj, iphd2, link_obj, ico_obj, &
            massp_obj, pxi_obj, peta_obj, pzta_obj, chip_obj, wate, chi, mv)
!
         write (*, *) 'celest3d: starting parcelc_obj'
!
         call parcelc_obj
      endif
!
!
      write (*, *) 'celest3d: starting vinit'
!
      call vinit
!
      write (*, *) 'celest3d: starting budget'
!
      call budget
!
!     ******************************************************************
!
!     generate an adaptive grid
!
!     *****************************************************************
!
      if (adaptg) then
!
         call gradf (nvtx, ijkvtx, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, &
            c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, &
            c3z, c4z, c5z, c6z, c7z, c8z, chic, gradxchi, gradychi, gradzchi)
!
         call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, &
            cdlt, gradxchi, gradychi, gradzchi, periodicx)
!
         call axisgrad (ibp1, jbp1, iwid, jwid, kwid, gradxchi, gradychi, &
            gradzchi)
!
!
!     set gradient at k=kbp2 equal to gradient at k=kbp1
!
         do j = 1, jbp2
            do i = 1, ibp2
!
               ijk = (kbp2 - 1)*kwid + (j - 1)*jwid + (i - 1)*iwid + 1
!
               gradxchi(ijk-kwid) = gradxchi(ijk)
               gradychi(ijk-kwid) = gradychi(ijk)
               gradzchi(ijk-kwid) = gradzchi(ijk)
!
            end do
         end do
!
         do n = 1, nvtx
!
            ijk = ijkvtx(n)
!
            chic(ijk) = sqrt(gradxchi(ijk)**2+gradychi(ijk)**2+gradzchi(ijk)**2&
               )/vvol(ijk) + 1.
!
         end do
!
         write (*, *) 'celest3d: calling smooth'
!
         nsm = 10
!
         call smooth (2, ibp2, 2, jbp2, 2, kbp2, iwid, jwid, kwid, nsm, chic, &
            gradxchi, periodicx)
!
         write (*, *) 'celest3d: starting meshgen'
!
         call meshgen (ibp2, jbp2, kbp2, x, y, z, wgrid, chic, vol, dt, taug, &
            numgrid, wate(1,7), wate(1,8), wate(1,9), wate(1,10), wate(1,11), &
            wate(1,12), wate(1,13), wate(1,14), wate(1,15), lama, lamb, lamc, &
            wate(1,1), wate(1,2), wate(1,3), delta, psi, tsix, tsiy, tsiz, etax&
            , etay, etaz, nux, nuy, nuz, wate(1,4), wate(1,5), wate(1,6), iwid&
            , jwid, kwid, ijktmp2, istep, iota, periodic, toroid, strait, rmaj&
            , rwall, btorus, q0, cdlt, sdlt, dz)
!
         write (*, *) 'celest3d: starting celindex'
!
         call celindex (2, ibp1, 2, jbp1, 2, kbp1, iwid, jwid, kwid, ijkcell, &
            ijkctmp, ncells, ijkvtx, nvtxkm, nvtx)
!
         write (*, *) 'celest3d: calling geom'
!
         call geom (ncells, ijkcell, nvtx, ijkvtx, iwid, jwid, kwid, ibar, jbar&
            , kbar, nvtxkm, cartesian, cdlt, sdlt, dz, x, y, z, c1x, c2x, c3x, &
            c4x, c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, &
            c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, vol, vvol, vvol_half, &
            ijktmp2, a11, a12, a13)
!
!
         write (*, *) 'celest3d: starting metric'
!
         call metric (ncells, ijkcell, iwid, jwid, kwid, x, y, z, tsix, tsiy, &
            tsiz, etax, etay, etaz, nux, nuy, nuz, sie)
!
!     relocate particles on the new mesh
!
         write (*, *) 'celest3d: starting parlocat_all'
!
         call parlocat_all
         write (*, *) 'celest3d: starting parlocat_all_obj'
!
!
         call parlocat_all_obj
!
         write (*, *) 'celest3d: starting parcelv_obj'
!
         call parcelv_obj (ncells, ijkcell, nvtx, ijkvtx, nrg_obj, iwid, jwid, &
            kwid, ijkctmp, itdim, iphead_obj, iphd2, link_obj, ico_obj, &
            massp_obj, pxi_obj, peta_obj, pzta_obj, chip_obj, wate, chi, mv)
!c
         call parcelc
!c
!     ADDED BY G. LAPENTA
         call parcelv (ncells, ijkcell, nvtx, ijkvtx, nsp, iwid, jwid, kwid, &
            ibar, jbar, kbar, dt, ijkctmp, anydrift, eandm, itdim, iphead, &
            iphd2, link, transpos, qom, ico, qpar, pxi, peta, pzta, up, vp, wp&
            , pmu, bxv, byv, bzv, wate, uptilde, vptilde, wptilde, mask, pixx, &
            pixy, pixz, piyy, piyz, pizz, qdnv, number, j12x, j12y, j12z, jmag, &
            jxtilde, jytilde, jztilde, jxs, jys, jzs,  &
            pixysp2, piyzsp2, pixzsp2, pixxsp2, piyysp2, pizzsp2, &
            pixysp1, piyzsp1, pixzsp1, pixxsp1, piyysp1, pizzsp1, uxsp1,uysp1,uzsp1,densp1,uxsp2,uysp2,uzsp2,densp2,0d0)
!
         call vinit
         call output
         call budget
!
      endif
!
!l    (write data file on disk)
!      call pltset (0)
!
!      call second (stim)
      write(*,*)'second-tlm',tlm
      call getjtl (itlm)
           write(*,*)'getjtl',itlm
      itlwd = 0
      t1 = stim
      t2 = t1
      trl = t1
      tmin = (0.02 + 0.005)*fibar*fjbar*fkbar
      tlm = itlm - tmin
      if (tlm <= 0) then
!        (time limit too small for problem)
         write (kpt, 9020) tmin
         write (kfm, 9020) tmin
9020 format('>>> error: time limit must exceed ',1p,e9.2,'seconds <<<')
         call exit (1)
!
!
!ll   3. time-step: integrate equations in time:
!        --------------------------------------
      endif
!

end subroutine initialization
